# THIS SCRIPT OUTPUTS THE REALLOCATION MAPS

#clear environment 
rm(list = ls())

library(readr)
library(tidyr)
library(haven) #to read in the Stata data file

require(tigris) #
require(tidycensus) #
require(sp)
require(rgdal) #
require(sf) #
require(raster)
require(stargazer) #
require(readxl)
require(haven)
require(zoo) #== for na.locf
require(snow) #
require(doSNOW) #
require(parallel)
require(foreach)
require(iterators)
require(foreign)
require(abind)
library(dplyr)
library(Rfast)

#set directories
WORK.DIRECTORY = "Z:/user/rih4/Solar Siting/JAERE Code Repository"
DATA.DIRECTORY = file.path(WORK.DIRECTORY,"Data JAERE")
WORK.OUT = file.path(DATA.DIRECTORY,"intermediate_output")
PLOT.OUT = file.path(DATA.DIRECTORY, "final_output")
dir.create(WORK.OUT);dir.create(PLOT.OUT, recursive=T)

reallocation.results <- read_csv(file.path(WORK.OUT,"Simulation_Results_all.csv")) %>% 
  mutate(zip = sprintf("%05d", zipcode)) %>%
  dplyr::select(-c(zipcode,NAME,STATE)) %>%
  relocate(zip)

#Reallocation maps
library(maptools)
US_zips <- readRDS(file.path(DATA.DIRECTORY, "FIPS Shapefiles/USA_Zip_Code_Boundaries_v104_zip_poly.rds"))
US_zips <- merge(US_zips, reallocation.results, by.x = "czip", by.y = "zip")

US_zips$original_systems[is.na(US_zips$original_systems)] <- 0
US_zips$interstate_env[is.na(US_zips$interstate_env)] <- 0
US_zips$interstate_env30[is.na(US_zips$interstate_env30)] <- 0
US_zips$intrastate_env[is.na(US_zips$intrastate_env)] <- 0
US_zips$intrastate_env30[is.na(US_zips$intrastate_env30)] <- 0
US_zips$interstate_env_energy[is.na(US_zips$interstate_env_energy)] <- 0
US_zips$interstate_env_energy30[is.na(US_zips$interstate_env_energy30)] <- 0
US_zips$intrastate_env_energy[is.na(US_zips$intrastate_env_energy)] <- 0
US_zips$intrastate_env_energy30[is.na(US_zips$intrastate_env_energy30)] <- 0


dots.original_systems <- dotsInPolys(US_zips, as.integer(US_zips$original_systems/10), compatible=FALSE)  # turn shapefile into dot density map (which is a SpatialPointsDataFrame)
dots.interstate_env <- dotsInPolys(US_zips, as.integer(US_zips$interstate_env/10), compatible=FALSE)  # turn shapefile into dot density map (which is a SpatialPointsDataFrame)
dots.interstate_env30 <- dotsInPolys(US_zips, as.integer(US_zips$interstate_env30/10), compatible=FALSE)  # turn shapefile into dot density map (which is a SpatialPointsDataFrame)
dots.intrastate_env <- dotsInPolys(US_zips, as.integer(US_zips$intrastate_env/10), compatible=FALSE)  # turn shapefile into dot density map (which is a SpatialPointsDataFrame)
dots.intrastate_env30 <- dotsInPolys(US_zips, as.integer(US_zips$intrastate_env30/10), compatible=FALSE)  # turn shapefile into dot density map (which is a SpatialPointsDataFrame)
dots.interstate_env_energy <- dotsInPolys(US_zips, as.integer(US_zips$interstate_env_energy/10), compatible=FALSE)  # turn shapefile into dot density map (which is a SpatialPointsDataFrame)
dots.interstate_env_energy30 <- dotsInPolys(US_zips, as.integer(US_zips$interstate_env_energy30/10), compatible=FALSE)  # turn shapefile into dot density map (which is a SpatialPointsDataFrame)
dots.intrastate_env_energy <- dotsInPolys(US_zips, as.integer(US_zips$intrastate_env_energy/10), compatible=FALSE)  # turn shapefile into dot density map (which is a SpatialPointsDataFrame)
dots.intrastate_env_energy30 <- dotsInPolys(US_zips, as.integer(US_zips$intrastate_env_energy30/10), compatible=FALSE)  # turn shapefile into dot density map (which is a SpatialPointsDataFrame)

dots.original_systems$simulation <- "Original"
dots.interstate_env$simulation <- "InterstateEnv"
dots.interstate_env30$simulation <- "Reallocation s.t. 30%"
dots.intrastate_env$simulation <- "Intrastate Env"
dots.intrastate_env30$simulation <- "Reallocation s.t. 30%"
dots.interstate_env_energy$simulation <- "Interstate Env Energy"
dots.interstate_env_energy30$simulation <- "Reallocation s.t. 30%"
dots.intrastate_env_energy$simulation <- "Intrastate Env Energy"
dots.intrastate_env_energy30$simulation <- "Reallocation s.t. 30%"


dots.original.interstate_env30 <- rbind(dots.original_systems, dots.interstate_env30)
proj_zip <- proj4string(US_zips)
proj4string(dots.original.interstate_env30) <- proj_zip

dots.original.intrastate_env30 <- rbind(dots.original_systems, dots.intrastate_env30)
proj4string(dots.original.intrastate_env30) <- proj_zip

dots.original.interstate_env_energy30 <- rbind(dots.original_systems, dots.interstate_env_energy30)
proj4string(dots.original.interstate_env_energy30) <- proj_zip

dots.original.intrastate_env_energy30 <- rbind(dots.original_systems, dots.intrastate_env_energy30)
proj4string(dots.original.intrastate_env_energy30) <- proj_zip


US_zips2 = st_as_sf(readRDS(file.path(DATA.DIRECTORY, "FIPS Shapefiles/USA_Zip_Code_Boundaries_v104_zip_poly.rds")))
US_zips2 = group_by(US_zips2, STATE)
US_zips2 = US_zips2 %>% st_buffer(0)
US_zips2 = US_zips2 %>% group_by(STATE) %>% summarize() 
US_states = as_Spatial(US_zips2)
library(rgeos)
US_states <- gSimplify(US_states, tol=0.01)
state.layer <- list("sp.polygons", US_states)
my.palette <- c("red", "blue")
point.size <- 0.4
dots.original.interstate_env30$simulation <- as.factor(dots.original.interstate_env30$simulation)
dots.original.intrastate_env30$simulation <- as.factor(dots.original.intrastate_env30$simulation)
dots.original.interstate_env_energy30$simulation <- as.factor(dots.original.interstate_env_energy30$simulation)
dots.original.intrastate_env_energy30$simulation <- as.factor(dots.original.intrastate_env_energy30$simulation)


map.interstate_env30 = spplot(dots.original.interstate_env30, "simulation", sp.layout = state.layer, col.regions = my.palette, cex = point.size, alpha=0.2, key.space = "right" , xlim=c(-124.74, -66.951), ylim=c(24.56, 49.39), par.settings = list(axis.line = list(col = 'transparent')))
map.intrastate_env30 = spplot(dots.original.intrastate_env30, "simulation", sp.layout = state.layer, col.regions = my.palette, cex = point.size, alpha=0.2, key.space = "right" , xlim=c(-124.74, -66.951), ylim=c(24.56, 49.39), par.settings = list(axis.line = list(col = 'transparent')))
map.interstate_env_energy30 = spplot(dots.original.interstate_env_energy30, "simulation", sp.layout = state.layer, col.regions = my.palette, cex = point.size, alpha=0.2, key.space = "right" , xlim=c(-124.74, -66.951), ylim=c(24.56, 49.39), par.settings = list(axis.line = list(col = 'transparent')))
map.intrastate_env_energy30 = spplot(dots.original.intrastate_env_energy30, "simulation", sp.layout = state.layer, col.regions = my.palette, cex = point.size, alpha=0.2, key.space = "right" , xlim=c(-124.74, -66.951), ylim=c(24.56, 49.39), par.settings = list(axis.line = list(col = 'transparent')))

proj4string(dots.original_systems) <- proj_zip
dots.original_systems$simulation <- "Installed"
dots.original_systems$simulation <- as.factor(dots.original_systems$simulation)
map.original_systems = spplot(dots.original_systems, "simulation", sp.layout = state.layer, col.regions = "red", cex = point.size, alpha=0.2, key.space = "right" , xlim=c(-124.74, -66.951), ylim=c(24.56, 49.39), par.settings = list(axis.line = list(col = 'transparent')))

#Save maps as PDFs
pdf(file.path(PLOT.OUT, 'Rplot_installed.pdf'), width=8, height=4)
map.original_systems
dev.off()
pdf(file.path(PLOT.OUT, 'Rplot_Interstate_Env30.pdf'), width=8, height=4)
map.interstate_env30
dev.off()
pdf(file.path(PLOT.OUT, 'Rplot_Intrastate_Env30.pdf'), width=8, height=4)
map.intrastate_env30
dev.off()
pdf(file.path(PLOT.OUT, 'Rplot_Interstate_Env_Energy30.pdf'), width=8, height=4)
map.interstate_env_energy30
dev.off()
pdf(file.path(PLOT.OUT, 'Rplot_Intrastate_Env_Energy30.pdf'), width=8, height=4)
map.intrastate_env_energy30 
dev.off()

